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CNJ [ Abstract. Current practice in stellar evolution is to employ one-dimen- 

£H ■ sional calculations that quantitatively apply only to a minority of the ob- 

^ | served stars (single non-rotating stars, or well detached binaries). Even in 

these systems, astrophysicists are dependent on approximations to handle 
^sO ! complex three-dimensional processes like convection. Understanding the 

structure of binary stars, like those that lead to the Type la supernovae 
used to measure the expansion of the universe, are grossly non-spherical 
and await a 3D treatment. 
I/"") ■ To approach very large problems like multi-dimensional modeling 

of stars, the Lawrence Livermore National Laboratory has invested in 
massively parallel computers and invested even more in developing the 
algorithms to utilize them on complex physics problems. We have lever- 
aged skills from across the lab to develop a 3D stellar evolution code, 
| Djehuty (after the Egyptian god for writing and calculation) that op- 

O ■ erates efficiently on platforms with thousands of nodes, with the best 

available physical data (opacities, EOS, etc.). Djehuty has incorporated 
all basic physics for modeling stars including an accurate equation of 
^ \ state, radiation transport by diffusion, thermonuclear reaction rates, and 

C$ ■ hydrodynamics, and we have begun testing it in a number of applications. 
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1. Why do stellar modeling in three dimensions? 

The best determinations of the size, age and composition of the universe are 
founded on measurements of stars, combined with a physical understanding 
based on one-dimensional (ID) computations. ID calculations quantitatively 
apply only to a minority of the observed stars (single non-rotating stars, or well 
detached binaries). Even in these systems, astrophysicists are dependent on ID 
approximations to complex three-dimensional (3D) processes like convection. 

Observations assure us that our best ID approximation of steady state 
convection is flawed, requiring ad hoc 'overshoot' corrections. Even worse, there 
is no real ability to model the time-dependent convection that may be critical 
in modeling helium flashes, Cepheid and Mira pulsations, or the nucleosynthesis 
in deep mixing events that follow thermal pulses. In addition, nearly half of the 
visible mass of the universe is incorporated into intrinsically 3D binary stars, 
many of which interact significantly at some stage of evolution. Understanding 
integral properties such as the chemical evolution, or the detailed pre-explosion 
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structure, of metric phenomena like type la supernovae awaits development of 
a true 3D stellar evolution code. 

Developing a 3D code capable of modeling stars is exceedingly challeng- 
ing. In spite of the inherent difficulty, the importance of obtaining a better 
understanding of stars has led a small number of groups to begin such work. 
Most current calculations are practically limited to of order 1 million zones, 
and cannot realistically handle whole stars (in 3D). Envelope convection has 
been simulated by modeling modest segments of a star (Clement 1993; Ludwig, 
Jordan, & Steffen 1994; Freytag, Ludwig, & Steffen 1996; Brummell, Hurlburt 
k Toomre 1996, 1998; Stein & Nordlund 1998; Skartlien, Stein, & Nordlund 
2000; Nordlund & Stein 2001; Porter & Woodward 2000). Quite reasonably, 
these simulations often lack physical processes pertinent to the core (nuclear 
energy production), or radiation transport, and the gravitational potential is an 
imposed condition. While these are important starts toward understanding con- 
vection in stars, it is clear that some problems (e.g. solar rotation) yield results 
that are dependent on the size of the segment that can be simulated (Robinson 
&; Chan 2001), necessitating whole star modeling. 

With the advent of massively parallel computers and the development of 
algorithms able to operate in such a partitioned environment, it became possible 
to consider developing a code, Djehuty, capable of following meshes with sizes of 
order 10 8 cells and with a more complete suite of physics. Djehuty is capable of 
modeling complete stars in three dimensions. It is based on the ALE (Arbitrary 
Lagrange-Eulerian) hydrodynamic method with a predictor-corrector Lagrange- 
Remap formalism, that is second-order accurate in both time and space. It 
supports energy transport, via a pair of diffusion equations, one for matter con- 
duction, and another for radiation. The mesh, constructed of multi-block logi- 
cally rectangular hexahedrons, is allowed to be non-orthogonal and can conform 
to the developing instabilities, providing a more accurate tracking. Djehuty per- 
mits domain decomposition for parallel operation, using MPI message passing. 

The physical data in Djehuty is identical to that used in the ID evolution 
code that provides initial structure models. It has an equation of state based 
on the work of Eggleton, Faulkner & Flannery (1973), as updated by Pols et 
al. (1995; hereinafter PTEH), that is very accurate for stellar compositions 
over a large range of masses (>O.5M ). Both Planck and Rosseland mean 
opacities for radiation diffusion come from OPAL (Rogers & Iglesias 1992) and 
are augmented by Alexander opacities (Alexander & Ferguson 1994) for the lower 
temperatures where molecules are significant. The nuclear reaction network 
includes reactions for hydrogen, helium, and carbon burning, and the nuclear 
reaction rates come from Caughlan & Fowler (1988). Nuclear energy production 
can also be computed by a tabular set of reactions. Neutrino losses are in tabular 
form from the work of Itoh et al. (1989, 1992, with errata). 

The gravity implementation is currently complete only for spherical stars, 
but is adequate to begin scientific investigation of a host of issues related to 
convection. In the coming year, we will implement a generalized gravity solver 
to study rapidly rotating or binary stars. We are still learning how to improve 
operational efficiency, but have begun a set of large scale, full physics, parallel 
runs aimed at our first science objective, the structure of convective cores in 
high mass main sequence stars. 
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2. The Physics in Djehuty 

2.1. An Astrophysical Equation of State 

As an analytic approximation, the PTEH equation of state (EOS) provides con- 
tinuous thermodynamic derivatives for hydrodynamic consistency. This EOS 
has been modified to reproduce tabulated OPAL values to an accuracy of better 
than 1% for the entire range of conditions expected in stars from 0.7 to 50 M Q 
(over their whole evolution). Stellar models as low as 0.5 M can be computed, 
with differences of only about 2% in their envelopes. 

A comparison of the difference between the OPAL tables and the PTEH 
EOS is shown in Figure 1. Over most of the temperature-density plane, the 
difference is under a percent. The upper left corner of the diagram is the region 
of complete degeneracy (relativistic and non-relativistic), important for white 
dwarfs, and the cores of red giants. While OPAL tables were not calculated 
there, the thermodynamics of this region is well represented by the Fermi rela- 
tions with which the PTEH EOS is in excellent agreement. The upper right hand 
corner is included in the OPAL tables but without the inclusion of relativistic 
effects. Here the PTEH code is presumably more accurate, and the difference 
shown as green and red contours arise from a limitation in OPAL calculations. 

Between the degenerate and non-degenerate regions there is a difficult region 
where Coulomb corrections to the pressure are very significant. This shows as 
a diagonal band of blue contours where the OPAL and PTEH EOS's differ by 
more than 1%. It is a region occupied by brown dwarfs, and very low mass stars 
( ^ 0.5 Mq), and we will improve agreement there if we decide to investigate such 
objects in the future. 

Successfully transferring a hydrostatic model to Djehuty (with minimal 
transients) requires that both codes utilize the same EOS. The PTEH EOS 
has been inserted into Djehuty and tested on a set of two and three dimensional 
problems as well as our 3D stellar model. In a simple hydrodynamics problem, 
the accurate EOS call can be a substantial component of the time. The time 
cost of using this accurate equation of state was minimized by implementing a 
sub-cycling procedure as well as a method improving the convergence rate of 
the algorithm. We have been able to optimize the run time of the EOS without 
sacrificing the accuracy of the results. From the first crude implementation of 
the accurate EOS, a five- fold speed up has been achieved. 

This EOS includes the formation of molecular hydrogen and ionization 
states of low-Z elements. While the EOS is able to operate on elements up 
to iron, it is implemented to consider only 1 H, 3 He, 4 He, 12 C, 14 N, 16 0, 24 Mg. 

2.2. Radiative Diffusion 

Energy transport via diffusion is an excellent approximation when the cells of 
the mesh are optically thick (t > 1 in a zone), a condition that is very well met 
throughout a star's interior. Djehuty uses a pair flux-limited diffusion equations 
to approximate energy transport. Included in the diffusion equation are sources 
and sinks that link the radiation equation to the hydrodynamic energy and 
momentum equations, as well as nuclear energy production, or neutrino losses. 
Djehuty operates in a two-temperature (2T) mode, in which the temperatures of 
the matter and of the radiation field are integrated separately, and the radiation 
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Figure 1. The thick contours mark regions of no difference between 
the analytic PTEH equation of state, and LLNL tabular values. The 
grey region marks regions where they differ by < 1%. The diagonal 
line tracks the temperature density profile of a solar model. 



is approximated as Planckian. Again, this is an excellent approximation within 
stellar interiors, and will enable us to transition to stellar atmospheres. These 
equations are: 

radiation : ^ = V. (— — V</>) + K p pc(aT 4 - 4>) + E$ (1) 
dT 

matter : P C V — = V. (D VT) + K p pc((p - aT 4 ) + E m . (2) 

Here a is the usual radiation constant, <f> is the radiation energy density per 
unit volume (which would be aT 4 in complete equilibrium), and E$,E m are 
contributions to the radiation and matter energies from compression, nuclear 
reactions etc. D is the heat transport coefficient by conduction alone. 

Within a star, the matter and radiation are so strongly coupled that the 
matter and radiation temperatures are essentially identical. This strong cou- 
pling is achieved because of the large value for the Planck opacity, k p . A table 
of Planck opacities has been constructed from data from the Opacity Project 
(Seaton et al. 1994) at high temperatures and of Alexander & Ferguson (1994) at 
low temperatures. Accurate Planck opacities will be important near the surface 
where the diffusion approximation becomes poor. 

The transfer of radiation is governed by the Rosseland mean opacity, K r . 
Currently Djehuty reads n r from a set of astrophysical opacity tables that were 
formed using OPAL (Livermore) opacities at high temperatures ( > 10, 000 K) 
and Alexander opacities at lower temperatures where molecules are important 
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(Alexander & Ferguson 1994; Rogers & Iglesias 1992). This is the same opacity 
set used by the ID evolution code. These astrophysical opacities will be the 
only option actually used in the stellar calculations done by Djehuty, but an 
analytic Kramer's form for the Rosseland mean opacity is available to facilitate 
comparisons with other codes. 

The radiation package of Djehuty has been given some simple tests to ensure 
accuracy. The first was a simple Marshak wave problem. A temperature source 
was applied to one end of a tube. With hydrodynamic motions turned off, and 
a constant opacity applied, the results can be compared directly with theory. 
Excellent agreement was found, providing direct evidence that the radiation 
diffusion package in Djehuty functions well. 

2.3. Thermonuclear Rates and Neutrino Losses 

Elements normally carried directly by the ID evolution code were selected to 
allow accurate tracking of principal energy generation reactions for hydrogen, 
helium and carbon burning over the whole range of stellar masses. This list of 
elements includes 1 H, 3 He, 4 He, 12 C, 14 N, 16 0, 24 Mg. These elements are passed 
to Djehuty as part of the mapping process, and form a standard 7-element set. 

The 7-element set is an excellent approximation for circumstances of sta- 
ble hydrogen burning in most stars. It includes the slower hydrogen burning 
reactions of the PP1 and PP2 chains as well as the CNO cycles. Over evolu- 
tionary time-scales, many reactions in the PP chain and CNO Cycles come to 
equilibrium, and do not need to be specifically included. 

The rate equations for the 7-element set are identical to those used in the ID 
code and Djehuty integrates those rates with the abundances to determine local 
energy production and transmutation rates. In this mode, energy production 
computed by Djehuty and the ID code are well matched. The rate equations 
used by Djehuty and the ID code include not just hydrogen burning, but the 
most important energy production rates from helium burning and carbon burn- 
ing stages. We have also included strong and weak Coulomb screening. 

Both Djehuty and the ID code have been written to permit easy expansion 
of the element set. Djehuty can also operate with an 18-element set suitable for 
situations where the hot CNO cycle operates. The program is now structured 
such that it is trivial to add new rates connecting existing elements. It is also 
relatively easy to increase the element sets when circumstances require. 

Operation of the new version is independent of the order in which the el- 
ements are specified in the input deck, and the solution of all rate equations 
is implicit. In tests, the hydrogen mass fraction can be changed by orders of 
magnitude in a single step. As part of these changes the EOS was modified to 
handle an arbitrary element set. 

In both element sets, the proton-proton chain is handled as follows: the 
proton capture on deuterium, p(p,/3 z^)D(p,7) 3 He, is assumed to be followed 



instantly by 3 He( 3 He,2p) 4 He and 3 He( 4 He,7) 7 Be(p, 4 He) 4 He. 

In the 7-element set, only the slower rates are included in the CNO cy- 
cle. The beta decay and proton capture on 13 C, 12 C(p,7) 13 N(/3,z/) 13 C(p,7) 13 N, 
are assumed to be instantaneous. The beta decay and proton capture on 15 N, 
14 N(p,7) 15 0(/3,^) 15 N(p,a) 12 C, are also assumed to be instantaneous, as are the 
beta decay and proton capture on 17 0, 16 0(p,7) 17 F(/?,z/) 17 0(p,a) 14 N. 
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All CNO rates are included in the 17-element suite making it suitable for 
the hot CNO cycle, including leakage into 19 F. They are: 

12 C(p, 7 ) 13 N, 13 N(/3,i/) 13 C, 13 C(p, 7 ) 14 N, 14 N(p, 7 ) 15 0, 15 0(/V) 15 N, 
15 N(p,a) 12 C, 15 N(p, 7 ) 16 0, 16 0(p, 7 ) 17 F, 17 F(/3, u) 17 0, 17 0(p,a) 14 N, 
17 0(p, 7 ) 18 F, 18 F(/?,i/) 18 0, 18 0(p,a) 15 N, and 18 0(p, 7 ) 19 F. 

The following helium burning reactions are also included: 4 He(2a, 7 ) 12 C, 
12 C(a, 7 ) 16 0, 14 N(a, 7 ) 18 0, and 18 0(a, 7 ) 22 Ne. In the 7-element set, the last 
reaction is assumed to happen instantaneously, and the mass fraction change is 
placed with all other heavy elements in 24 Mg. 

Finally, the following reactions are included for beginning advanced stages 
of massive star evolution: 12 C( 12 C, 7 ) 24 Mg and 16 0( 16 0, 7 ) 32 S. 

As a result, the 7-element set uses 12 reaction rates and the 18-element set 
evaluates 23 rates. These rates are all evaluated from equations following the 
form of Harris et al. (1983) and returned as reactions per mole per second. How 
the reactions are used is defined by three numbers returned for each rate that 
specify the number of input particles, of output particles, and the number of 
particles in reaction. This permits the code to include appropriate factors for 
true multi-body rates like the triple-alpha rate. It also allows two-body rates, 
like 12 C(p, 7 ) 13 N, to be correctly evaluated while reducing the number of protons 
by 2 in the approximation where the reactions 13 N(/3,i/) 13 C(p, 7 ) 14 N are assumed 
instantaneous. 

Each rate also includes the net energy produced from the rate (less any 
neutrino energy). This information is used along with the rates to derive the 
energy production rate in ergs per gram per second. A table for neutrino losses 
(Itoh et al. 1989, 1992) from Compton, bremsstrahlung, and pair production 
processes has been included, and is initialized along with the basic data for the 
EOS routine. The temperature and density are then used to interpolate in this 
table to provide a neutrino loss rate. 

Finally, Djehuty can be quickly transformed by providing the rates in tab- 
ular form. We have generated a table using rates associated with the 7-element 
set and shown good correspondence in energy production and transmutation 
rates. The tabular method can be implemented without recompiling Djehuty, 
and lends itself to quick tests. In the future we anticipate using it to simulate 
chemical rates for modeling giant planet atmospheres. 

2.4. Gravity 

At present we have implemented a simple algorithm to compute the gravity of 
a spherical object. Between hydrodynamic cycles, Djehuty integrates over all 
cells of the mesh to determine a mass-radius relation. In parallel operation, this 
requires sharing a relatively small amount of information among the processors 
to construct the relation for the entire star. 

Once the spherical mass distribution has been found, the gravitational ac- 
celeration can be simply calculated at each node. This method is limited to 
spherical objects, but its implementation identified much of the process that 
must occur for a more general treatment to be implemented. 

As an extension of this self-consistent spherical gravity, it is possible to add 
a point mass outside of the problem, and study tidal modification of the model. 
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This approximation is appropriate only when the star is sufficiently centrally 
condensed that the portions of the star that are distorted are of modest mass. 

2.5. Hydrodynamics 

Djehuty uses an ALE (Arbitrary Lagrange-Eulerian) hydrodynamic method, in 
which an explicit hydrodynamic step is followed by an advection step. In stable 
regions, or regions of large-scale coherent motion where the Lagrangian mesh 
is modestly deformed, no advection is necessary, and the code is essentially 
Lagrangian. In regions where shear develops, it is necessary relax the mesh, and 
permit material to move between cells. 

To promote hydrodynamic accuracy, and time centering for other physical 
processes, the hydrodynamic step is split in half. The first sub-step is a simple 
explicit step. This is followed by a second sub-step in which information from 
both the original position, and information from the half (sub) step position, are 
used in a predictor-corrector formalism. This allows a time-centered inclusion 
of radiation transport and energy production with a hydrodynamics step that 
is second-order accurate in both time and space. In Djehuty, vector quantities 
like position and velocity are zone-centered, while scalar quantities (density, 
temperature ...) are cell-centered. By splitting the hydrodynamic step, the 
node and cell quantities have the same time. 

3. Model Generation 
3.1. The Initial Models 

We developed a ID stellar evolution code for use as a platform to test physics, 
and to provide structure information to Djehuty for constructing 3D models. 
Djehuty reads the ID stellar evolution models directly, and converts the model 
into its own variable set. The evolution code follows the work of Eggleton (1972) 
and has a fully implicit algorithm that includes mesh motion during a time step. 
This speeds the calculation and improves its accuracy, but results in an unusual 
centering of the thermodynamic variables. To minimize any transients, software 
was also written to adjust (and relax) the variables to values appropriate for 
cell-centered equations used by Djehuty. 

As an example the ID code is capable of following a 4 M & model from a 
pre-main-sequence configuration through main sequence, giant, and thermally 
pulsing stages to becoming a white dwarf (Figure 2). At any stage of evolution, 
a restart file can be written that is usable by either Djehuty or the ID code 
itself. 

As a first test of Djehuty, a pre-main-sequence model of a star approaching 
the main sequence was tested. At this point, nuclear reactions are negligible 
(E nuc / E grav < 0.005), the composition homogeneous, and the star is almost en- 
tirely stable. This model was selected so as to test the accuracy of our scaling 
and mapping algorithms as well as to test various 3D mesh structures. 

A second model that was selected for an initial test case was of a star that 
has just reached the main sequence. It has a convective core of about 0.4 R & 
(2.8 x 10 10 cm), driven by the strong temperature sensitivity of CNO nuclear 
burning. In this model (Figure 3), carbon has been converted into nitrogen, 
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Figure 2. An HR diagram shows the temperature luminosity evo- 
lution of a 4Mq star from the pre-main-sequence to a point where 
its remnant core has become a white dwarf. The model positions are 
shown by open circles. 

the core 3 He has been destroyed, and the central hydrogen mass fraction has 
fallen from 0.7 to 0.676. This model has been used to test the nuclear energy 
generation rates as implemented in Djehuty. 

3.2. Mesh Structure 

In order to reduce the problem of tiny zones (and time steps), inherent in a mesh 
using spherical coordinates, we are utilizing a technique in which the sphere is 
constructed from 7 logical blocks. There is a central block, and 6 adjacent blocks, 
one attached to each face of the central block. The mesh in each of the outer 
blocks is morphed to form a spherical tile. The central block can remain as a 
cube or be morphologically distorted from near rectangular cells at the center to 
a sphere at its surface. Preliminary numerical experiments suggest that allowing 
the central block to retain most of its cubic shape results in a better-behaved 
structure, and smaller deviations from hydrostatic equilibrium (particularly at 
the corner nodes of the central cube). 

As an example (Figure 4), we show a sample logical structure in which 
the central block (0) is constructed from a mesh with 50x50x50 cells (51x51x51 
nodes). There are then 6 additional blocks (1-6), having 50x50x100 cells. In 
each case, the axis oriented in the radial direction has 100 zones (101 nodes). 
Blocks 1, 2, and 3, attach to the faces on which the logical variables i, j, and k 
are 1, while blocks 4, 5, and 6 attach to the faces where i, j, and k have their 
maximum values. 
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Figure 3. The composition structure of an early main sequence star 
of4M„. 



Logical structure of Djehuty Mesh 



A 1.6 million zone mesh 
sphere is constructed of 7 
logical blocks. The central 
block contains 50X50X50 
cells, and each face is 
connected to a block of size 
100X50X50. 
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Figure 4. Logical structure of Djehuty mesh. 



10 



Bazan et al. 




Figure 5. The structure of a representative mesh with 590K zones 



The logical (surface) blocks that attach to the central cube must (at present) 
have the same mesh size as the face to which they adjoin, but the radial compo- 
nent may have any length. These surface blocks are transitioned from a cubic 
shape to spherical segments (Figure 5). An initial mesh is generated with a 
standard size. It can then be tailored to a particular stellar model (radius and 
radial zoning structure). The scaling must not just match the radius of the star, 
but have a radial spacing to resolve important gradients. The radial structure 
of the ID mesh is used as an initial guide for scaling the 3D mesh, but the user 
has an ability to bias the mesh through a parameter. 

This additional flexibility is provided by producing a scale factor that is 
parameterized as 'grid_scale_alpha'. The scale factor, S r , was just defined as 

S r = r a (3) 

where r is the fractional radius of a node in the mesh. Setting grid_scale_alpha 
to a value of 1 results in a direct mapping of the radial structure from the ID 
model to the 3D mesh. A value of grid^cale_alpha greater than 1 concentrates 
the mesh in the core, and makes the envelope sparser. Alternatively a value less 
than 1 enlarges the core and refines the envelope (Figure 6). 

Once the mesh has been scaled, the ID model is mapped into it. The cells of 
the small meshes that we have used for our initial testing are of such a size that 3 
or 4 of the ID radial zones contribute to the properties of the cells. Furthermore, 
the cells near the core will not (in general) have a radial orientation. A simple 
mapping of zone-centered values led to a cumulative error in total stellar mass 
that was unacceptably large. This mapping problem was solved by subdividing 
our 3D cells and integrating the ID model across the smaller volume elements 
(Figure 7). The segment masses and volumes were then summed to determine 
the average cell density. Other properties were then mass-averaged among the 
segments. 
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Figure 6. In a coarse 3D mesh, 3 to 4 radial zones from the ID model 
may contribute to the properties of a cell. 
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Figure 7. Comparison of ID model to Djehuty mesh. 



4. Model Settling 



Figure 8 gives an example of the initial velocity perturbations that result from 
the mapping of a 4M model to a 3D mesh. The model is generated with a 
fixed boundary condition, and run for 5 cycles to an elapsed time of 0.074416 
seconds. The grey scale ranges from to 500 cm/s. 
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Figure 8. Mesh structures. The initial velocity transients are mini- 
mized when the transition region is smoother. 
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5. Conclusion 

Djehuty is now an operational code, and papers are in preparation on stellar 
simulations that have been made with it (see Figure 3 of Eggleton et al. in these 
proceedings) . Highlights of these calculations include the first 3D simulation of 
a type II supernova. The nickel produced in the silicon-burning region of these 
stars is an important energy source in the observed light curves. In SN1987a, 
it appeared much sooner that expected from a ID model. Our 3D calculation 
shows the development of an instability causing tendrils of nickel to extrude 
into surrounding material, where it appears early, as seen by observers. We 
have also done a first simulation of a nova. Novae are explosive events that 
occur in binary-star systems where hydrogen is being torn from a normal star 
and deposited on to a white dwarf companion. In our first 3D calculation, the 
hydrogen ignition propagated from a number of ignition points, and created a 
shock that accelerated the surface to near 3000 km/s. When the shock emerges 
at the surface, it is heated to temperatures in excess of 1 Kev, creating an X-ray 
pulse. In this first calculation, the white dwarf mass is a bit smaller that those 
thought appropriate for most novae, and the hydrogen accretion was symmetric 
over the surface of the star. A more realistic hydrogen distribution in which 
the material is settled on to the star as a high velocity stream will be a future 
calculation. Simulating stable stars, like our 4M© main sequence model, has 
been remarkably challenging, but with new boundary conditions the numerical 
stability problems that precludes long runs seem to have been resolved. 
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